Jet trails and Mach cones: The interaction of microquasars with 

the ISM 



Yoon, D. 1 , Morsony, B. 1 , Heinz, S. 1 , Wiersema, K. 2 , Fender, R.P. 3 , Russell, D.M. 4 , & 

Sunyaev, R. 5 ' 6 

ABSTRACT 

A sub-set of microquasars exhibit high peculiar velocity with respect to the 
local standard of rest due to the kicks they receive when being born in supernovae. 
The interaction between the radio plasma released by microquasar jets from such 
high- velocity binaries with the ISM must lead to the production of trails and bow 
shocks similar to what is observed in narrow-angle tailed radio galaxies and pulsar 
wind nebulae. We present a set of numerical simulations of this interaction that 
illuminate the long term dynamical evolution and the observational properties of 
these microquasar bow shock nebulae and trails. We find that this interaction 
always produces a structure that consists of a bow shock, a trailing neck, and 
an expanding bubble. Using our simulations to model emission, we predict that 
the shock surrounding the bubble and the neck should be visible in H Q emission, 
the interior of the bubble should be visible in synchrotron radio emission, and 
only the bow shock is likely to be detectable in X-ray emission. We construct an 
analytic model for the evolution of the neck and bubble shape and compare this 
model with observations of X-ray binary SAX J1712. 6-3739. 
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Introduction 



It is now well established that accretion onto black holes can lead to the production 
of powerful jets, both in the case of AGN, in which case we call the object a radio galaxy, 
and in the case of X-ray binaries (XRBs), in which case we call the object a microquasar. 
We have since learned that even neutron star X-ray binaries as well as some accreting white 



dwar f s behave in the sam e way (IFender et al.l 120041 ; iMigliari fc Fender! l2006t iTudose et al. 
20091 ; iKording et alJliooi ). 



When powerful jets run into their environment, they produce bright external shocks 
(called hot spots or working surfaces) and generally inflate pockets of relativistic plasma 
that emit diffuse synchrotron emission. If the black hole is stationary, these pockets take on 
roughly ellipsoidal shapes and are called radio lobes or cocoons. They are surrounded by the 
interstellar or intergalactic gas that was occupying the volume prior to inflation, compressed 
into a narrow shell. This picture is most easily observed in the case of X-ray cavities in 
galaxy clusters inflated by AGN jets ( Rafferty et al. 20061. and references therein), but has 
also been discovered around a few micro quasars (jGallo et al.ll2005t lHao &: Zhang! 120091 ). 



However, when the black hole is moving at appreciable speed with respect to its en- 
vironment, the ram pressure of the headwind can dramatically alter the dynamics of the 
outflowing radio plasma, sweeping it back and generating a bow shock a head of the movin g 
black hole. This has long been known in the case of radio galaxies (e.g. iMiley et al.lll972l ). 



which, depending on the angle at which the plasma is bent back, are called narrow or wide 
angle tail sources (generically, these sources are also called "bent doubles"). 

Given that the black holes at the core of microquasars are born in supernova explosions, 
and give n that supernovae c an impose significant kick velocities on the compact objects they 
produce, iHeinz et al.l (120081 ) suggested that a similar phenomenon to "bent doubles" should 
occur around a sub-population of microquasars: the high-velocity tail of the population of 
low-mass X-ray binaries (LMXBs). It is now known that several mic roquasars are movin g: 
through the ISM at relative velocities in excess of v ext > 100 km s -1 (IMirabel et al.ll200ll ). 



The phenomenology of these trailed microquasars should be broadly similar to their 
supermassive AGN equivalent, leading to the production of a low-surface brightness trail of 
synchrotron emitting relativistic plasma, and a brighter bow shock nebula. While the mode 
of inflation is fundamentally different, the end product should also be very similar to pulsar 
wind bow-shock nebulae that are formed by interaction between strong wind from pulsars 
with significant kick velocities in the interstellar medium. The general bow shock fe atures 
have been well observed in H„ emission for PSR B1957+20 (IStappers et al.l 120031 ). PSR 



B0740-28 f lStappers et al.l 120021 ). PSR J0437-4715 (Bell et al.l Il995h . and PSR B2224+65 
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(the "guitar" nebula) (IChatteriee fc Cordesll2002h . Also, G359. 23-0.82 which is powe red by 
PSR J 1747-2958 shows bright head and X-ray & radio "trails" JGaensler et al.ll2004h . 



The first candidate of a bow shock around such a trailed source, the LMXB SAX 



J1712.6-3739 , was dis covered in an H Q image by IWiersema et al.l (120091 originally found b y 
in 't Zand et al.lll999l ). It is broadly consistent with the predictions of iHeinz et al.l ( 120081 ). 
Future searches for large scale nebula and diffuse synchrotron emission around known LMXBs 
are needed to further test the predicted existence of this population of sources. 



The model developed in IHeinz et al.l (120081 ) was purely analytic and a number of im- 
portant aspects of the evolution of such trailed microquasars was left open. Numerical 
simulations are needed to investigate the product i on of these sources in more detail and to 
confirm the qualitative predictions of IHeinz et al.l ( 120081 ). 



In this paper, we model the evolution of the XRBs moving through the ISM by using 
3-dimensional hydrodynamic simulations and to derive a more detailed dynamical model for 
the large scale evolution of these sources. 

A key aim of this study is to develop quantitative diagnostics that can be used to derive 
important constraints on the core parameters of microquasars from observations of trailed 
microquasars, such as their age, their relative velocity with respect to the local standard 
of rest (uxrb), and the jet power from observational parameters like the opening angle of 
the bow shocks, the size of the hot bubble at the terminus of the trail, and the brightness 
of the shock, trail, and bubble. In addition, our simulation results can be used to predict 
the brightness of H a , bremsstrahlung, X-ray, and radio synchrotron emission and to design 
targeted observational searches for these sources. 

This paper is organized as follows: In §2, we present the numerical method and the initial 
conditions for our numerical study. In §3, we discuss the evolution of XRBs, scaling relations, 
and observational expectations. In §4, we compare our results with new observations of the 
LMXB SAX J1712. 6-3739. Finally, in §5 we summarize our results. 



2. Technical Description 



2.1. The code 



Simulations are carried out in 3 dimensions with the FLASH 2.4 hydrodynamic code 
( iFryxell et al.l |2000| ). which is a modular, adaptive mesh refinement code. It solves the 
Riemann problem using the piecewise-parabolic method. The code is formally accurate to 
second order. The gas is modeled with either an adiabatic equation of state with index 
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7=5/3 or with radiative line cooling. Radiative cooling from an optically thin plasma is 
implemented in FLASH adopting a piecewise-power law approximation that evaluates a 
plausible fit to A(T) within the temperature range of 4.4 x 10 3 < T < 10 8 K 



2.2. The jet nozzle 

In order to simulate the injection of collimated, superso nic jets into the gr id, we employ 



a numerical "nozzle", as first developed and described in ( IHeinz et al.ll2006l ): An internal 
inflow boundary of cylindrical shape placed at the location of the XRB, injecting fluid with 
prescribed energy, mass, and momentum flux to match the parameters we chose for the jet. 

For reasons of numerical stability, we impose a slow lateral outflow with low mass flux 
in order to avoid complete evacuation of zones immediately adjacent to the nozzle due to 
the large velocity divergence at the nozzle. The injection of energy and mass due to this 
correction is negligible. 

We generally follow the prescri ption for jet injection used in previous simulations of AGN 



jets described in (IHeinz et al.ll2006l ). but keep the location of the XRB fixed in space, instead 
letting the external medium stream by at velocity — v ext . Unlike our AGN simulations, we 
do not impose a random jitter on the jet axis in this case. 

We chose to inject the jet at an internal Mach number of 10. While our adiabatic 
simulations are strictly scale free, we picked a set of fiducial dimensions for our box that 
lead to the following natural scaling: For computational feasibility, we chose a jet velocity of 
v jct = 3 x 10 9 cms _1 . The jet is turned on initially and continues to inject material for the 
entire length of the simulation. 

The simulations were carried out with Adaptive Mesh Refinement in order to capture 
the large dynamic range required, ensuring that the nozzle is resolved with at least 8 cells 
across. For our fiducial scaling, the maximum resolution for the standard model is about 
0.156 pc near the jet nozzle. 



Based on the estimated power of the jet in Cyg X-l from I G alio et al.l (120051 ). the jet 



power in our simulation is set to a constant value of W- Jet = 10 37 ergs s" 1 in our fiducial 
scaling. Note that since the Cyg X-l is known to be one of the most powerful XRB sources, 
sustaining the hard X-ray state for about 90% of the time. The jets in our model might be 
more powerful than those of typical LMXBs. However, with the exception of our simulations 
with cooling, our simulations are scale free and can thus be adjusted easily to other parameter 
combinations. 
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2.3. Initial conditions 



We placed the XRB in a moving medium inside a box large enough that boundary 
conditions never influence the dynamics. We varied the velocity of the gas relative to the 
XRB to be f xrb = 30 km s _1 , 100 km s _1 and 300 km s _1 (see Tabled]). The lowest velocity 
case represents a typical LMXB, given that the LMXB velocity dispersion is 37kms _1 , while 
the largest represents the most extreme case plausible (a marginally unbound source). The 
simulations we ran for this study are listed in Tab. [TJ 



The fiducial ISM pressure we use is Pq = 3 x 10 -12 ergs cm -3 following ICoxl ( 120051 ) . 
and we use an ISM number density of n ISM = 1 cm -3 , giving a sound speed in the ISM of 
c s = 17.3 km s _1 . 

Typically, simulations were carried out to 1 Myr in our fiducial scaling (much longer 
than the dynamical evolution of the jet, and long enough for the quasi-steady state of the 
bow shock and the self-similar solution we will discuss below to be established). 

Table [1] also includes one case with our standard parameters and radiative cooling (with 
vxrb = 300 km s _1 ), the hvc case. We used this run to verify that the morphology of our 
simulations are not strongly affected by cooling and to properly model emission (see §3.3). 
A further 3 simulations with cooling, fitl, fit 2 and fit 3, were carried out to attempt to match 
E a observations of SAX J1712.6-3739 (see §4.1) 

Because of the self-similar nature of the flow, the simulations lend themselves to a 
staggered refinement scheme. As we will describe below, the radio plasma forms of a large 
scale bubble and a neck connecting the XRB to this bubble. We can therefore de-refine the 
simulation on cylinders on increasing radius around the axis of propagation of the XRB and 
again de-refine with increasing distance from the XRB along the axis of propagation. We 
decrease the refinement by a factor of two for every power of two increase in radius and 
distance. 



Table 1. Parameter of the simulations 



Name 


^XRB 


Maximum Resolution 


Luminosity 


Radiative Cooling 


lv 


30 km s" 1 


0.312 pc 


10 37 erg/s 


no 


hit 


100 km s" 1 


0.156 pc 


10 37 erg/s 


no 


hv 


300 km s" 1 


0.156 pc 


10 37 erg/s 


no 


hvc 


300 km s- 1 


0.156 pc 


10 37 erg/s 


yes 


fitl 


100 km s" 1 


0.039 pc 


8.47 X 10 33 erg/s 


yes 


fit2 


200 km s- 1 


0.039 pc 


8.47 X 10 33 erg/s 


yes 


fit3 


300 km s- 1 


0.039 pc 


2.29 X 10 35 erg/s 


yes 
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Results 



3.1. Morphological Evolution 



Heinz et al.l (120081 ) presented a simple, analytic model for the long term dynamical 
evolution of trailed microquasars. Our simulations allow us to move beyond the initial 
heuristic model and describe the dynamics of microquasar driven bow shock nebulae in 
detail. 



The initial evolution of the system follows the prediction of iHeinz et al.l (120081 ) almost 
exactly: The jets inflate a single large, roughly spherica l bubble that expan ds following the 
self- similar solution for a continuously driven bubble by ICastor et al.l ( 119751 ). surrounded by 
a thin shell of swept-up ISM. 

Initially, the expansion velocity of the bubble is much faster than the space velocity 
of the XRB, vxrb, implying that the XRB remains roughly at the bubble's center. As the 
expansion slows down, the XRB begins moving towards the shell and eventually breaks out 
of the expanding bubble. 

At this point, the dynamic pressure of the ISM due to the XRB's velocity causes the 
jets to bend backwards, and a trail of radio plasma is created behind the XRB, connecting 
it to the bubble. The radio plasma released by the XRB continues to inflate the bubble, 
akin to a balloon inflated by a straw. The bubble, therefore, continues to expand spherically 
even though the binary can be far outside the bubble. 

As it propag ates, the XRB is driving a bow shock into the ISM. As predicted in 



Heinz et al.l (120081 ). dynamical instabilities develop at the contact discontinuity between the 
bow shocked ISM and the backflow along the radio trail, dissipating some of the backflow 
energy along the channel. 

The dynamical evolution described above is shown in Fig. [T]as a time sequence of density 
slices through the cent er of the simulation box. R elativistic, axisymmetric simulations of a 
pulsar wind nebula in iBernstein &: Hughes! ( 120091 ) created a similar structure with a bow 
shock connected to an expanding bubble. This indicates that, away from the jet or wind 
source, the evolution of XRBs and pulsar nebulae are nearly identical, except for the scale. 



3.2. Scaling relations 



Based on the observed evolution of the trail and bubble, we can construct a simple 
analytic model for the inflation of bow shock nebulae by microquasars (and, by extension, 
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Fig. 1. — Density slice through our 300 km s _1 simulation at 10,000 years (top), 100,000 
years (middle) and 1 Myr (bottom). The XRB is located toward the left of the images. 
At 10,000 years (top) the XRB is surrounded by a spherical bubble, although the XRB is 
nearing the left edge. At 100,000 years (middle), the XRB has broken out of the bubble and 
created a distinctive shape. The jets, directed up and down, are curved back at the bow 
shock (curved shock at left) and directed back to the large bubble to the right, where they 
terminate and power the bubble expansion. A neck connects the binary and bubble, and is 
surrounded by an oblique shock. At 1 Myr (bottom), the shape is similar to that at 100,000 
years, except that the neck is more elongated, with a smaller opening angle. 
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pulsar wind bow shock nebulae) that can be tested directly against the simulations. This 
will allow us to formulate analytic expressions for observables. 

The jets initia lly inflate a nearly spherical cavity that evolves as a wind driven bubble 
(ICastor et al.lll975l ). The bubble expands as 

, / T \ 1/5 



where L is the source luminosity, po is the ISM density and C\ is a constant that equals 



25 
14tt 



for an adiabatic index of T 
this is 



5/3. In terms of the fiducial parameters for our simulations, 



R b (t) = 13.1 parsec ( — ^ t 3 J 5 (2) 
\PohJ 

where L 37 = lo37 ^ gs -i , Pm = L67xl0 -Z g cm -a , and t 5 = 100[00 f 0yoars . The expansion velocity 
of the bubble is 



t*(f) = fC^ 5 f V5 t^/ 5 = 76.7 km s" 1 f ^ V5 1" 2/5 (3) 
\PoJ \PohJ 

and the pressure inside the bubble is 

P b (t) = P {\M 2 b - |) (4) 

where Mj, = , = is the Mach number of the spherical shock and Po is the ISM pressure. 

V7-Po/po 

Assuming at least a moderately strong shock, the pressure goes to 

P b (t) « \P,Ml « ^C x 2/5 L 2 / 5 ^ /5 r 4 / 5 « 8.25 x lO" 11 erg cm- 3 L^ 5 po^ 5 " 4/5 (5) 

Initially, the expansion of this bubble is faster than the velocity of the XRB. However, 
as the expansion of the bubble slows while the velocity of the source remains constant, the 
source eventually breaks out of the bubble at a time tbrcak- The breakout time scales with 
source velocity v s as 

t break = c^(-\'\:^ 



\Po 



U 7 \ 1 ' 2 ( V X \-5/2 



185, 000 years ( — ] ( T ) ' (6) 

\PohJ VlOOkms" 1 / 
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The radius of the bubble at the breakout time therefore scales as 

/r\l/5 / r \ 1/2 / r \ 1/2 _ 3/2 

/ \ .3/5 ^1/2 / ^ \ _s/9 ~ / -^37 \ / \ d/Z 



100 km s 



-i 



(7) 



After breakout, the source continues to power bubble expansion. A bow shock develops 
in front of the source which bends the jets back in the opposite direction of the source motion. 
A low-density channel remains which connects the source to the bubble. Jet material flows 
back through this channel, adding energy to the interior of the expanding bubble. So long 
as the volume of this channel remains small compared to the volume of the bubble, the 
expansion rate of the bubble continues to be described by eqn. [fl the self-similar equation 
of a continuously powered bubble. 

As seen in Fig. [TJ the resulting shape is an expanding spherical bubble connected to 
the source by a thin neck. The neck consists of a shock surrounding a narrow cavity filled 
with jet material. At any point x along the neck, the shape is described by w±, the width 
of the inner cavity, and u>2, the distance from the mid plane to the outer edge of the shock. 
In the frame of the binary, the cavity has a fixed shape and is in pressure balance with 
the surrounding shocked material. The flow of material can therefore be described using an 
adiabatic equation of state, the Bernoulli equation and mass continuity: 

P x = apl (8) 

V + 1 £k = 

2 (7 - 1) P. 



PxV x wl = C (10) 

where v x , p x and P x are the velocity, density and pressure of material in the cavity at position 
x, 7, the adiabatic index, is 5/3, and a, b and c are constants. In term of a, and with 7 = 5/3, 
b = a 3 / 5 (p ^) 2/5 and c = L a" 3 / 5 ^ 2 )" 2 / 5 . 

Rearranging eqns. [8] to [10] we solve for W\ and find that 
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Making the approximation P x <C p v^ (i.e., the source velocity is supersonic), this reduces 
to 



Wl = 5 -l/4 2 l/ 2jL l/ 2a 3/10 (poU 2 ) -3/10p-3/10 = C2 p-mo (12 ) 

The pressure will be set by the jump conditions for an oblique shock, with 9 and /3 the angle 
with respect to the direction of motion of the inner and outer edges of the shock, respectively: 



tan(0) 



M 2 sin^ 



B-l 



tan/3(l + ^M 2 



Po 



2jM< 



sin 



- M 2 sin 2 0) 
2 /3-( 7 -l) 



(7 + 1) 



(13) 
(14) 



where Mo = v s / a/7-Po/Po is the Mach number of the binary relative to the ISM, tan# = 
dwi/dx and tan (3 = dw2/dx. Making the approximation that 9 and j3 are small, this reduces 
to 



Substituting we find 



A(fM 2 /3 2 - |) (16) 



/ 1 (tci/ga)- 1 °/ 3 + i 



5M 2 ' ^ 




dun = 3(C 2 10/3 - Ppwf 3 ) / 20C 2 10/3 
~ dx~ M P < /3 \Powf 3 

This differential equation can then be integrated numerically to find wi(x) and this in 
turn can be used to find W2(x). 

As the pressure drops along the neck, however, it eventually reaches a minimum value 
equal to the pressure in the expanding bubble. Beyond this point, the angle of the shock 
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ia constant such that the post-shock pressure is equal to the P b , the pressure in the bubble, 
which is set by the expansion rate of the bubble v The incoming velocity towards the shock 
will be v s — Vb, because the bubble is expanding towards the source, giving a minimum shock 
angle of 



= ^- = ^* (19) 
v s — v b dx 



The angle of the inner surface is then approximately 



M|/3 2 - 1 _ d Wl 
±M$/3 ~ ~dx 



$ = T„o = -T 1 (20) 



where M 2 = (v s — v 5) / a/ 7P0 / Po is the Mach number of the shock. In our simulations, this 
asymptotic limit dominates and the shape described by eqn. [TS] only describes a small region 
near the jet source. Note that the width of the neck where it reaches the minimum value of (3 
is described by eqn. [12] with P x = Pb, giving wi(Pb) = C<iP b 3 ^ 10 - This equation depends on 
the normalization constant a in the adiabatic equation (eqn. which in practice depends 
on the amount of mixing between jet and ISM material. The width at this point also depends 
on the radius of curvature of the bow shock, which is not taken into account in this analytic 
model Therefore, we use a as a free parameter to get the proper fit for the width of the neck 
where the minimum value of /3 is reached. 

Figure [2] shows density slices of the 300 km s _1 , 100 km s _1 and 30 km s _1 simulations 
(models hv, mv and lv) at 1 Myr. The white lines are results of our analytic model for the 
bubble size, 101(2;) and 102(2;). In the first two cases, the predicted shape of the inner and 
outer edge of the neck are a good fit to the simulations. In the 30 km s _1 case, the XRB is 
still inside the spherical bubble at 1 Myr. 

The outer shock angle (3 decreases with time as the pressure in the bubble and the 
expansion velocity drop. The width of the neck where it meets the bubble is approximately 



W 2 ,meet = P(v s t - R b ) + W 1 (P b ) (21) 

w 2 , meet = {v s t - \ Vb t) + C 2 (|P M 6 2 ) - 3/1 ° (22) 

^b 

Asymptotically, this width goes to W2 t meet = 0v s t = v b t = vRb- Therefore, the ratio of 
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the bubble radius to the neck width approaches a constant ratio of 3/5, and the width of 
the neck is always smaller than the radius of the bubble. 

The volume of the neck scales asymptotically as 



K cck = > s t (fifc) 2 = \*Cl'\ (£) 2 V/s = 3.23 x 10- cm3 {^- T ) t?" 

(23) 

while the volume of the bubble scales as 



V bub = \kRI = fvrCf (±) 3/5 t^ = 2.75 x 10 59 cm 3 tf (24) 

the ratio of the volumes is therefore 



^ = — ^ = ^\ ( ^V'V = 0.117 ( " s A ( ^V V5 tf (25) 
V bub 100 R b 20 1 s \poJ VlOOkms-V \ Pm ) 5 1 ' 

The volume of the neck and bubble will eventually become equal when the length of the 
neck (v s t) is about 11.1 times the bubble radius. The width of the shock where it meets the 
bubble is about W2 )me et = §-R&, so the angle when the volumes are equal is = 27/500 = 0.054 
However, the angle of the shock around the neck cannot be less that /3 ~ 1/M , so the volumes 
become equal before the shock becomes weak only if Mo > 18.5. Our assumption that the 
volume of the neck is small will hold until bubble expansion starts to become marginally 
sonic, unless the source has a very high Mach number relative to its surroundings. 

While our simulations have not run long enough to probe the sub-sonic regime of bubble 
expansion, it is worth speculating about the late state evolution of trails and bubbles. Given 
that the bounding pressure of the bubble will be dominated by the internal pressure of the 
ISM, the expansion velocity of the bubble will drop below the self-similar value for an energy 
driven bubble (with R oc t 1 / 3 ). In addition, the ISM will no longer be strongly compressed 
into a narrow shell. Since the neck and bow shock will maintain their stationary shape, one 
should expect that the late state evolution of a bow-shock/trail nebula will eventually lose 
the terminating bubble and the trail pressure will eventually approach th e ISM pressure, 



consistent with the late state structure of the trail proposed in iHeinz et al.l (120081 ) . 
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Fig. 2. — Density slice of 300 km s" 1 (upper), 100 km s" 1 (lower left) and 30 km s _1 
(lower right) simulations at 1 Myr. The white curves trace the analytic solution for W\ and 
u> 2 , the neck shape, and the circle is the size of the bubble from eqn. [TJ Values of a used to 
fit the bow shock width are 10 31 for the 300 km s _1 simulations and 10 33 for the 100 km 
s™ 1 simulation. In the 30 km s _1 simulation the XRB (marked with an x) is still inside the 
bubble, but significantly offset from the center. 
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3.3. Observable properties of microquasar bow shock nebulae and trails 

3.3.1. H a emission 

For XRBs moving supersonically through cold or warm ISM, the compressed, hot gas in 
the bow shock will produce collisionally excited line emissi on, most importantly (as well 



as the classic spectrum of nebular lines like [OIII]). In fact, IWiersema et al.l ( 120091 ) detected 



the nebula of SAX J1712. 6-3739 in H Q . Our models support the jet-ISM interaction scenario 



( IHeinz et al.ll2008l ) and the numerical simulations show an apparent quantitative agreement 



with the observed results for H a emission. 



To calculate the H„ emission, we first determine the ionization balance in the shocked 



gas, using the MAPPINGS III code ([Sutherland fc Dopital Il993l ); the code uses a time 



dependent algorithm for accurate equilibrium balance calculation. The typical temperature 
of the shocked shell is estimated to be around 10 5 — 10 6 K in the adiabatic case, hence 
the gas is inferred to be nearly fully ionized at the shocked region. Figure |3] (bottom left 
panel) shows the surface brightness in H a for the model hv. Note that all simulations can 
be scaled with XRB jet luminosity by reducing the predicted surface brightness, nebula size, 

^37) 



and simulation time by a factor of {L%-^ 1 ' 2 



The bubble is relatively bright in H a . The shock around the nec k is also visible in Hr, i n 



fig. [3j A similar feature has been identified with SAX J1712. 6-3739 (IWiersema et al.ll2009l ). 



We discuss the morphological similarity between that source and our simulations in §4.1. 

Radiative cooling has little influence on the dynamical evolution of the XRBs. The 
upper panels of fig. |3] reveal that overall morphologies for two models (with and without 
radiative cooling) are broadly similar. Radiative cooling leads to a significantly thinner shell 
of shocked material around the bubble and neck, giving the appearance of a slightly narrower 
neck. Because of the cooling-induced contraction of the gas, the shell also appears to develop 
some irregularities in shape, though the overall shape of the hot (radio) plasma inside the 
trail and bubble occupies essentially the same volume. 

However, the radiative cooling has a pivotal role for H Q , because of the strong temper- 
ature sensitivity of the ionization balance and thus the recombination line emission. This is 
because the cooling time scale is comparable to the dynamical time scale of XRBs, leading 
to significantly lower temperatures in the shell than in a purely adiabatic simulations. 

The bow shock immediately surrounding the XRB is very dim because the temperature 
is too high for strong H Q emission. As cooling becomes important down stream, the temper- 
ature at the shock quickly drops to the range of a few 10,000 K. As a result, the emission 
becomes stronger by about a factor of 50. 
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Because H a emission depends non-trivially on temperature, it is not straight forward 
to express brightness predictions semi-analytically. Figure [3] is therefore specific to our 
simulation using fiducial parameters only. 
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Fig. 3. — The density contour map for the model without the radiative cooling (upper left) 
and with the radiative cooling (upper right) at 100,000 yr (v xr b s = 300 km s _1 ). The bottom 
panels shows H Q emission for each case (log scale). The predicted surface brightness can be 
scaled to different XRB jet luminosities by reducing the surface brightness, size scale and 
simulation time by a factor of L^ 2 . 
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3.3.2. Radio emission 

Given that the trail itself should be filled by magnetized, relativistic plasma released by 
the jets, they will emit synchrotron radiation. As the plasma travels along the jet, particles 
will cool both adiabatically and due to synchrotron losses. As is well known from studies of 
AGN jets and radio lobes, this introduces a cutoff to the electron energy distribution and to 
the synchrotron spectrum. 



Following, e.g., iHeinz fc Begelmanl (119971 ). the cutoff frequency, as a function of travel 



time along the trail, will be given by 



p 



i/i 

(26) 



I¥tUb{^ dt 
with an associated cutoff frequency of 

max Anm e c* ' 1 ' 

Given that the flow through the trail assumes a quasi-steady state behind the XRB, 
we can simply integrate this equation through single frames to lowest order to derive the 
cooling frequency as a function of position along the trail. The result is shown in fig. HI For 
our fiducial parameters, the cooling frequency is estimated to fall within 10 13 — 10 15 Hz, well 
above the radio band. 

This implies that synchrotron emission should be a good tracer of these trails at all 



freque ncies, not just at low frequencies. This is in contrast to the estimates in IHeinz et al. 



( 120081 ). who suggested that cooling could be important along the trail. The reason for 



this differenc e is the sign i ficant backflow velocity along the trail, which was left as a free 



parameter in IHeinz et al.l (120081 ) . and which allows radio emitting particles to traverse the 



trail without significant losses. 

Figure shows synchrotron surface brightness for the 300 km s _1 simulation at 100,000 
years. Surface brightness for synchrotron emission is strongest at the front-edge where the 
jets are bent by ram pressure. The terminal radio bubble is also bright and emits the bulk 
of the total radio luminosity. 

In fig. O the typical brightness temperature inside the bubble is estimated to be of order 
a few K in model hvc at 0.1 Myr at a frequency of 1 GHz. Since the synchrotron emission 
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Fig. 4. — Synchrotron cooling frequency for X-ray binaries at 100 km s~ l (left panel) and 300 
km s -1 (right panel) at 1 Myr. Top images indicate the density map with logarithmic scale. 
The vertical magenta lines indicate edges of the bubble area (marked with the solid circle). 
The cutoff frequency is calculated by integrating along the two blue lines and assuming a 
steady state. 



depends only on the pressure, the typical values of brightness temperature in other models 
(^xrb=30, 100 km s" 1 ) are essentially the same. 

If the intensity of synchrotron emission is high enough, synchrotron self-absorption will 
occur and the emission will drop out at low frequencies, proportional to u 5 ^ 2 regardless of 
the electron power index. However, in our particular simulation, the optical depth at 1 GHz 
is quite transparent (about r = 10~ 2 ) so this effect would be negligible. 

The bright temperature can be extrapolated into a few 10~ 3 K in submm wave bands 
that might be detectable by a new generation detectors such as SCUBA-2 or ALMA. In the 
case of the XRBs in our galaxy, assuming that the distance would be a few kpc, the angular 
size of the bubble is large enough to be resolved by these detectors. 

Since the pressure inside the bubble decreases with time, the surface brightness will also 
decrease with time and synchrotron radiation will be easier to observe at an earlier XRB age 
if the bubble is resolved. 

The synchrotron emission from the bubble can easily be predicted analytically. For a 
power-law distribution of electrons, N(j)dj = A^~ v d^\ if we assume the power-law index 
(p) is 2.5, the total emissivity per unit volume per unit frequency can be calculated as, 
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Fig. 5. — Surface brightness of radio synchrotron emission for the frequency of 1 GHz at 
100,000yr {v X rb = 300 km s" 1 ). 
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Fig. 6. — Surface brightness of radio synchrotron emission for the frequency of 1 GHz at 
30,000yr {v XRB = 300 km s" 1 ). 
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tsync — C 3 P 15 / 8 , (28) 

where P is the pressure in the bubble and C3 is approximately 1.44 xlO -17 and depends on 
the ratio of electron to magnetic pressure. Throughout this discussion we will assume that 
the plasma is in equipartition (with equal energy density in electrons and magnetic field), 
making the estimated fluxes upper limits. 

The pressure of the hot bubble can be estimated from the shock jump conditions. If we 
set the adiabatic index to 5/3, the pressure jump condition in the shell of the bubble are 
given by eqn. |5j 

Figure [7] shows that the inferred post-shock pressure is consistent with the simulation 
results, although the data from the simulation have a slightly smaller values. The small 
discrepancy is negligible in estimating the surface brightness. 

The surface brightness can be derived from the emissivity (eqn. 128]) with post-shock 
pressure and the path length of the line of sight in bubble. Hence the maximum value of 
path length is Rb u b(t), and the surface brightness becomes 

J ~ C 3 ( 27 r 2/5 j 2/5 3/5 4 /5 \ 15/8 1/5 J l/5 -1/5.3/5 

isync ~ ^— I 1 ^° ) 1 ^° 

/ f \ -9/10 

« 1.37 x 10- 3 mJy arcsec" 2 x (L 37 ) 19/2 ° {p m f 1/m (^J (29) 

This theoretical calculation is in good agreement with numerical data shown in fig. [HJ The 
relation confirms that the synchrotron emission decreases over time, implying that younger 
XRBs are more easily detected in the radio. 

XRB bubbles in other galaxies will be difficult to resolve with radio observations, but 
may be detectable as unresolved sources, at locations different from the X-ray point source. 
If we set the distance to the XRBs to be D, then we can estimate the unresolved flux as, 



Flux sync « ^(^-Cl /6 L 2 "pTt' A A C^L^p-^D- 2 
y 16vr V100 1 Ho J 1 Ho 

" - 04mJ >' x '« 2,,20 (^ 2,/ ii^) 3/IO (i^) 



3/10 

(30) 
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This flux includes only the bubble area, and neglects the synchrotron emission from near the 
jets. Because we do not resolve the base of the inner regions of the jet in our simulations, 
and because microquasars are generally highly variable, we cannot make any quantitative 
statements about the relative flux from the inner jet and the bubble. Note that the total 
flux from the trail and bubble increases with time, indicating that older XRBs will be easier 
to detect than young ones in unresolved observations. 



O simulation (mv) 
Scaling relation 





Fig. 7. — The solid line indicates that post shock pressure and density derived from scaling 
relation and diamond marks the results from simulation. 

The state-of-the-art radio observatories are well suited for this type of objects; the 
eVLA has been upgraded to have the excellent sensitivity, frequency coverage, and imaging 
capability that allows for good quality radio spectra resolved over the source, and the WSRT 
APERTIF that has very large field of view at 1.4 GHz, allowing for deep continuum surveys 
capable of detecting trailed nebula. 
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Fig. 8. — The solid line indicates the theoretically expected surface brightness in synchrotron 
emission and diamond indicates the results from simulation. 

3.3.3. Free-free emission 

For the expected temperature range in the bow shock, we can expect the gas to be 
mostly or completely ionized, which will give rise to the emission of Bremsstrahlung. Figure 
IH] shows a simulated radio image in free- free emission for the case of wxrb =300 km s" 1 . 
The emission is edge-brightened, given that it comes exclusively from the compressed ISM 
shell. 

The temperature at the shell of the bubble can be calculated by the post-shock density 
and pressure as 

bub ~ k Pl 400 Gl k Po 1 ' 

where /i is mean molecular mass. If we assume fully ionized gas inside the shell, the free-free 
surface brightness through the shock can be approximately written as 



I Jf « 3.91 x 10" 20 ergs s-W^-W 1 X {p m f (r^\ 
where we assume a strong shock with p s hock ~ 4p for simplicity. 



(31) 
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Free-free emission is dimmer than the synchrotron emission in the bubble, and is ex- 
pected to be undetectable or marginally detectable at ^=lGhz. 

For the shock temperatures in our simulations (10 4 — 10 6 K), optical emission is domi- 
nated by cooling lines rather than free-free continuum emission. However, for a fast source 
the temperature in the bow shock may initially exceed 10 6 K, in which case free-free emission 
would be the dominant coolant. 




Fig. 9. — Surface brightness of free-free emission for the frequency of 1GHz at 100,000 yr 
(yxrbs = 300 km s _1 ). 



3.3.4. X-ray 

For large enough space velocities, the temperatures at the stagnation point of the bow 
shock can be sufficiently high to produce X-rav s. Fig. [TOl shows a simu lated Chandra X-ray 



observation produced using the XIM program jHeinz fc BruggerJl2009h for the 300 km s 



simulation at a distance of 5 kpc with a 100 ks exposure at 10 6 years. 

For the fiducial parameters, there is a significant flux only from the bow shock, which 
has a temperature of about 2 x 10 6 K. However, it may be difficult to distinguish this flux 
from the emission from the XRB. There is a small amount of X-ray emission from the neck, 
but it is only marginally distinguishable from the background, and fades as you go farther 
down the neck. 
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At early times, the expanding bubble produces a detectable X-ray flux. Figure [TT] shows 
a simulated Chandra X-ray observation for the 300 km s -1 simulations (hvc) at a distance 
of 8 kpc with a 100 ks exposure at 30,000 years. Although it has a low surface brightness, 
the bubble and neck can be seen, in addition to the bright bow shock, in a smoothed image. 

A shock will only be detectable in X-rays if its temperature is above about 5 x 10 5 
K, which corresponds to a shock velocity of about 100 km s _1 for the parameters in our 
simulations. Assuming a strong shock, the shock temperature is 



T s ~ 6.79 x 10 3 ^ xM s 2 xf ( ^ -) (32) 

s \3 x io~ 12 y Vi.67 x io- 24 y V ; 

The bubble expansion falls below about 100 km s" 1 at about 36,000 years. After this, 
the bow shock is still visible, if the XRB is moving fast enough, but the bubble and neck are 
unlikely to be detected. The bow shock is always visible, but only for XRBs moving faster 
than about 100 km s" 1 . 

The synchrotron emission from relativistic particles does not produces a detectable X- 
ray flux. 



4. Comparison with SAX J1712. 6-3739 



While the numerical results and the analytic approximations derived from them are 
primarily predictive, we can, at the very least, directly apply them through comparisons to 
the one known XRB with a bow shock nebula. 

SA X J1712.6-3739 is so fa r the only XRB found to display a prominent H a bow shock 
nebula (IWiersema et al.ll2009i ). Following the original discovery of the bow show, a deep 
VLT observation of the source was obtained to confirm the detection and look for further 
structure in the H-alpha map. This observation is shown in Fig. | 



4.1. Data reduction 

The data were taken using the Focal Reducer and low dispersion spectrograph (FORS2) 
on the Very Large Telescope (VLT, Chile), using the H Q +83 interference filter. As part of a 
larger programme (ESO programme ID 385.D-0100, PI Russell), we obtained 34 exposures of 
60 seconds each on 16 April 2009, under excellent seeing conditions (average 0.7 arcseconds). 



2 4 6 8 

x [pc] 




Fig. 10. — Simulated Chandra X-ray image (top) and surface brightness (bottom) of the 
bow shock for the 300 km s _1 simulation at 1 Myr. 100 ks exposure at 5 kpc and flux is 
integrated from 0.3 to 3 keV. The maximum flux in the bow shock is about 1 count/pixel. 
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Fig. 11. — Simulated Chandra X-ray image (top) and surface brightness (bottom) of the 
bubble for the 300 km s _1 simulation (hvc) at 30,000 years. 100 ks exposure at 8 kpc and 
flux is integrated from 0.3 to 3 keV. 
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The data were reduced using standard procedures in IRAF and combined. Details on 
the reduction and flux calibration of these data will be reported in a forthcoming publication 
(Russell et al. in prep.). The FORS2 data cover a much larger sky area than the EFOSC2 
data reported in Wiersema et al. (2009) at a better resolution, and have significantly better 
signal to noise, allowing us to search for fainter features and both small and large scale 
structure. 



4.2. Model comparison 

The n ew VLT observat i on co nfirms the presence of a bright bow shock nebula, as first 



reported in I Wiersema et al.l (120091 ) . In addition, the image suggests the presence of a roughly 
circular shell of H a emission with a radius of 76" which connects to the previously observed 
linear H Q features near the XRB. While the level of patchy background H a emission makes 
a firm identification difficult, we interpret this cavity or shell as the bubble inflated by the 
backflow through the channel. The estimated location of the bubble and neck are sketched 
in the lower panel of figure [12j The distance from the end of the neck to the center of the 
bubble is about 190" . 

We estimate the half opening angle of the bow shock seen in the H Q emission to be 
about 15°, which can be converted into the projected oblique shock angle of /3 = 0.261799. 
From the observed ratio of bubble size to neck length, our analytic model would predict a 
neck opening angle of about 19°, somewhat larger than the observed value. Assuming the 
bubble is expanding supersonically, we can set a lower limit on the Mach number of the 
source from eqn. [T9] of 



M s > ^Y~M b > 4.8 (33) 

For our assumed sound speed of 17.3 km s _1 , this implies the velocity of the source is 
v s > 83 km s" 1 . The source could also have a component of motion towards or away from 
us which is not constrained by the opening angle of the neck. 

Given the distance of the object of D ~ 7kpc and the approximate angular size of the 
putative shell of 76", we estimate the physical radius of the bubble to be -Rbubbie ~ 2.6 pc. 
In order to check the consistency between the observation and our model, we carried out a 
simulation with v s = 100 km s _1 . Using the scaling relation of bubble radius from eqn. (j2J) 
and a standard ISM density of 1 cm -3 , the appropriate luminosity is L = 8.5 x 10 33 ergs s" 1 



and the age of the bubble would be approximately 70,000 years: 




) 



1/3 



Rb 



) 



5/3 



t S Ax ~ 70, 000 yr 



2.57pc 



(34) 



where L 



L 



8.5e33 8.5X10 33 ergs s" 1 ' 



In terms of source power and ISM density, the velocity of SAX J1712. 6-3739 can be 



A higher energy (or lower ambient density) implies a higher velocity required of the source 
by eqn. [35] to have the same neck length at a given bubble size. If the motion of the binary 
is not perpendicular to our line of sight, but rather inclined at an angle a, then the true 
opening angle of the source is (3 true = (3sma and the velocity of SAX J1712. 6-3739 would 
have to increase by 1/ sin a. Eqn. [35] is therefore a lower limit on the source velocity for a 
given set of parameters. Based on the velocity analysis, a proper motion of SAX J1712.6- 
3739 can be estimated as a few mas/yr. This proper motion scale will be observable from 
the LSST. If the velocity is measured directly, it would help to constrain our models and 
allow us to determine the jet power in terms of just the ambient density or vice versa. 

In figure [T2], the Hq, emission in the neck is brighter than that in the expanding bubble 
while we expect it to be relatively dim from our fiducial simulations. 

It is possible that this is due to the inclination angle for the moving XRBs with respect 
to the line of sight. If the source is moving at a substantial angle with respect to the plane 
of the sky, then the amount of neck material we are looking through increases while the 
thickness of the spherical bubble stays the same. As a result, with some large inclination 
angle, the emission of the neck can be brighter than that of the bubble. 

It is also possible that radiative cooling in the shell around the bubble has moved much 
of the gas to temperatures too low to emit in H Q , while the neck is still at a sufficiently 
high temperature to emit, given that the bow shock and Mach cone have a higher Mach 
number and pressure than the bubble and shell. Because radiative cooling breaks the scale 
invariance of our simulations, a more detailed investigation would require the construction 
of a large grid of simulations at different source powers and ISM temperatures, which would 
be well beyond the scope of this paper. 

It is also possible that the relative brightness is simply a result of the complicated 
morphology of background emission surrounding the observed H a nebula. The surrounding 



estimates to be 




1/3 



(35) 




Fig. 12. — SAX J1712. 6-3739 image from the VLT FORS2 H a data (upper panel), and same 
image with location of the neck and bubble highlighted (lower panel). 
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emission is clearly non-uniform. A more detailed and quantitative investigation of the new 
VLT observations is forthcoming in a separate paper. 

The difficulty in interpreting the complex structure of the observation notwithstanding, 
we carried out three simulations with the initial conditions expected to produce a size and 
shape close to SAX J1712. 6-3739 to test whether we could match the observed shape and 
relative brightness of the neck and bubble. The parameters of these simulations are listed in 
table 1 as fitl, fit 2 and fit 3. In figure [T3"| the top panel shows the results when the ejected 
energy is 8.47 x 10 33 ergs s _1 and the velocity of the XRB is 100 km s _1 (model fitl). There 
is bight H Q emission in the bow shock and dimmer but still bright emission from the bubble. 
However, the emission in most of the neck is relatively dim compared to the shell. 

In order to inspect the inclination angle effect, we ran simulation fit 2 with a higher 
velocity of 200 km s" 1 and display it with an inclination angle of 60° (middle panel). In 
this case the emission is brightest in the bow shock, part of the neck and the leading edge 
of the bubble, with somewhat dimmer emission from the trailing edge of the bubble. This 
is somewhat closer to the bright neck seen in SAX J1712. 6-3739, although the location of 
bright neck emission does not match the observation exactly. The width of the neck relative 
to the bubble size is about the same as the 100 km s _1 case, and there is still a large gap 
in the neck with very little emission. 

Finally, we ran simulation fit 3 with a larger power of L — 2.29 x 10 35 ergs s _1 , and higher 
velocity velocity of ^xrb = 300 km s" 1 (bottom panel). For the given bubble size, the age is 
reduced to 23,000 years (eqn. |3"4"1) . Note that because the emission is significantly brighter in 
this model, the intensity scale for the bottom panel is increased by a factor of 10 compared 
to the other two panels to avoid color saturation. In this case, there is no emission from 
the bow shock, due to the high temperature caused by the faster source. There is bright 
emission in part of the neck, where the shock has become cool enough to emit in H a , and 
there is about equally bright emission from the bubble edges. By varying the source velocity, 
it is possible to move the location of bright neck emission from the bow shock (for model 
fitl) far down stream in the neck (for model fit 3) and it could be moved even farther away 
from the XRB for a faster source velocity. 

These three simulations differ substantially in their relative surface brightness distribu- 
tions of H Q emission in bow shock, neck, and shell, showing that it should be possible to 
construct a reasonable fit to the observations with sufficient fine tuning. 
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Fig. 13. — Hq, images with different initial conditions; L = 8.47 x 10 33 ergs s _1 , t>xRB 
1 00 km s" 1 and edge-on at 69,000 yr (top panel), L = 8.47 x 10 33 ergs s -1 , v X rb = 200 km s" 1 
and inclination angle of 60° at 69,000 yr (middle panel), L = 2.29 x 10 35 ergs s _1 , ^xrb = 
300 km s _1 and edge-on at 23,000 yr (bottom panel). Note that the intensity scale of the 
bottom panel is 10 times higher than the other two panels. 
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5. Conclusion &c Summary 

By way of numerical simulations, we study the morphological evolution and dynamics 
of microquasars with high space velocity, the functional equivalent of narrow-angle tail radio 
galaxies for accreting stellar mass black holes and neutron stars. Such dynamic X-ray bina- 
ries initially create a powered, spherically expanding shock/bubble, the expansion of which 
eventually decelerates to less than the velocity of the XRB relative to the ISM. 

When the source reaches the edge of the bubble, it establishes a well defined structure 
with a bow shock around the source at one end, a spherical bubble at the other, and a neck 
connecting the bow shock and bubble. The shock angle around neck, /3, decreases as the 
binary moves away from the bubble, but it cannot be less than 1/Mo, where Mo is the Mach 
number of the binary relative to the surrounding medium. 

The shocks from all three components creates H a emission that should be detectable in 
narrow band imaging. At early times, less than about 30,000 years in our simulations, X-ray 
emission will be detectable around the bubble and neck, but this rapidly fades as the shock 
temperature decreases. At later times, X-ray emission is only detectable from the bow shock 
and only for fast sources (v s > 100 km s" 1 ). 

Synchrotron radio emission from the bubble is bright and covers a large area, but emis- 
sion from the neck is significantly weaker. There is also strong synchrotron emission from 
near the binary, but it only covers a small area and thus would be difficult to detect without 
high resolution. Bubbles created by galactic XRBs should be visible in surveys of diffuse 
continuum emission. Binaries in other galaxies may also be detectable with a resolution such 
that the beam size is approximately the size of the bubble. The free-free emission is dimmer 
than the synchrotron emission and is generally not detectable. 

New H Q observations of SAX J1712. 6-3739 presented here (see §H]) show both a strong 
bow shock and tentative evidence for a spherical shell trailing the XRB, consistent with the 
neck and bubble morphology predicted by our simulations. We are also able to constrain the 
Mach number of the XRB relative to the background ISM to be M s > 4.8, corresponding to 
a velocity of v s > 83 km s _1 in the plane of the sky for a sound speed of 17.3 km s _1 . 

Extending our results into pulsar bow shock nebulae (PWBN) is plausible, since bow 
shock structures of PWBN are well describe d by the simulations presented in this paper. Fo r 



example, the H a images of PSR B2224+65 flChatteriee fc Corded 12002 : ICordes et al.lll993j l 



called "Guitar nebula" for its peculiar shape, shows well developed structures of trailing neck 
and spherical bubble. We will present models specific to PWBN in a future paper. 
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